if !d.name eq 'PS' then begin
    device,xsize=18,ysize=12,yoffset=3
    !p.thick=2 & !x.thick=1 & !y.thick=1
end

@eu_setup
@rhoprof

ns=93D
ng=50D
rsinth1=fltarr(m,l)
d_fr_dphi=fltarr(n,m,l,ns)
d_fth_dphi=fltarr(n,m,l,ns)
dsin_fphi_dth=fltarr(n,m,l,ns)
d_fr_dth=fltarr(n,m,l,ns)
dr_fth_dr=fltarr(n,m,l,ns)
dr_fphi_dr=fltarr(n,m,l,ns)
;
fphis=fltarr(n,m,l,ns)
fths=fltarr(n,m,l,ns)
frs=fltarr(n,m,l,ns)
;
; curl components
cfr=fltarr(n,m,l,ns)
cfth=fltarr(n,m,l,ns)
cfphi=fltarr(n,m,l,ns)
;
; total fields
fphi=bx
fth=by
fr=bz
;
; large-scale fields
fphim=total(fphi,1)/double(n)
fthm =total(fth,1)/double(n)
frm  =total(fr,1)/double(n)
;
; small scale fields
for i=0, n-1 do begin
 fphis[i,*,*,*]=fphi[i,*,*,*]-fphim[*,*,*]
 fths[i,*,*,*]=fth[i,*,*,*]-fthm[*,*,*]
 frs[i,*,*,*]=fr[i,*,*,*]-frm[*,*,*]
endfor


;SMALL SCALE HELICITY COMPONENT

for it=0,ns-1 do begin
print,'Computing total helicity',it
 for k=0, l-1 do begin
  for i=0,n-1 do begin
   dsin_fphi_dth[i,*,k,it]=xder_curl(reform(smooth(sin(theta)*fphis[i,*,k,it],3)),theta)       
   d_fr_dth[i,*,k,it]=xder_curl(reform(smooth(frs[i,*,k,it],3)),theta)
  endfor
 endfor

 for j=0,m-1 do begin
  for i=0,n-1 do begin
   dr_fth_dr[i,j,*,it]=xder_curl(smooth(r[*]*fths[i,j,*,it],3),r)   
   dr_fphi_dr[i,j,*,it]=xder_curl(smooth(r[*]*fphis[i,j,*,it],3),r)
  endfor
 endfor
;

 for k=0, l-1 do begin
  for j=1,m-1 do begin
   rsinth1[j,k]=1./(r[k]*sin(theta[j]))
  endfor
 endfor
 rsinth1[0,*]=0.
 rsinth1[m-1,*]=0.


 for k=0, l-1 do begin
  for j=1,m-1 do begin
   for i=0,n-1 do begin
     cfr[i,j,k,it]=rsinth1[j,k]*dsin_fphi_dth[i,j,k,it]
     cfth[i,j,k,it]=-dr_fphi_dr[i,j,k,it]/r[k]
     cfphi[i,j,k,it]=dr_fth_dr[i,j,k,it]/r[k]-d_fr_dth[i,j,k,it]/r[k]
   endfor
  endfor
 endfor

endfor

cht=total((cfr*frs + cfth*fths + cfphi*fphis),1)/double(n)
chtm=total(cht[*,*,ns-ng:ns-1],3)/double(ng)
mhel=fltarr(m,l)

for k=0, l-1 do begin
 mhel[*,k]=chtm[*,k]/(4e-7*rh[k])
endfor

save,filename='mhelicity.sav',r,th,x,y,mhel
end
